module nh_case6_mod
  ! Gravity waves
  use const_mod

  implicit none

  private

  public nh_case6

  real(8), parameter :: Rd      = 287.04
  real(8), parameter :: Rv      = 461.5
  real(8), parameter :: Cp      = 1004.5
  real(8), parameter :: Cv      = 717.5
  real(8), parameter :: omega   = 0!2 * pi / 86164.
  real(8), parameter :: p0      = 100000
  real(8), parameter :: ps      = p0
  real(8), parameter :: p_top   = 0
  real(8), parameter :: t0      = 300
  real(8), parameter :: u0      = 0!40
  real(8), parameter :: gamma   = 0.005
  real(8), parameter :: dpt     = 10
  real(8), parameter :: Lz      = 20000
  real(8), parameter :: radius  = 6371229
  real(8), parameter :: gravity = 9.80616d0
  real(8), parameter :: lon_c   = pi
  real(8), parameter :: lat_c   = 0

contains

  subroutine nh_case6(lon, lat, z, rho, u, v, w, pt, mr, zs, rc, g, omg, radi, zt)
    real(8), intent(in )           :: lon
    real(8), intent(in )           :: lat
    real(8), intent(in ), optional :: z
    real(8), intent(out), optional :: rho
    real(8), intent(out), optional :: u
    real(8), intent(out), optional :: v
    real(8), intent(out), optional :: w
    real(8), intent(out), optional :: pt
    real(8), intent(out), optional :: mr
    real(8), intent(out), optional :: zs
    real(8), intent(out), optional :: rc ! Rayleigh damping coefficient
    real(8), intent(out), optional :: g
    real(8), intent(out), optional :: omg
    real(8), intent(out), optional :: radi
    real(8), intent(out), optional :: zt
    
    real(8) :: t, p, f, S, N, r, RR, hs, theta, kappa
    
    N  = 1.e-2
    RR = radius / 3.
    kappa = Rd / Cp
    
    r  = spherical_distance(lat_c,lon_c,lat,lon,radius)
    hs = merge( 0.5 * ( 1. + cos( pi * r / RR ) ), 0._8, r<RR )
    
    if( ( present(rho) .or. present(u) .or. present(v) .or. present(w) .or. present(pt) .or. present(mr) .or. present(rc) ) &
    .and. .not. present(z) )then
      stop 'Need z in nh_case1 while acquiring rho, u, v, w, pt, mr, rc'
    endif
    
    if( present(z) )then
      S     = gravity**2 / ( Cp * N**2 )
      p     = p0 * ( ( 1. - S / T0 ) + S / T0 * exp( -N**2 * z / gravity ) )**( Cp / Rd )
      theta = T0 * exp( N**2 * z / gravity ) + dpt * hs * sin( 2. * pi * z / Lz )
      t     = theta * ( p / p0 )**kappa
    endif
    
    if (present(radi  )) radi  = radius
    if (present(rho   )) rho   = p / ( Rd * t )
    if (present(u     )) u     = u0 * cos(lat)
    if (present(v     )) v     = 0
    if (present(w     )) w     = 0
    if (present(pt    )) pt    = theta
    if (present(mr    )) mr    = 0
    if (present(zs    )) zs    = hs
    if (present(rc    )) rc    = 0
    if (present(g     )) g     = gravity
    if (present(omg   )) omg   = omega
    if (present(zt    )) zt    = Lz / 2.

  end subroutine nh_case6
  
  ! spherical distance on unit sphere
  function spherical_distance(lat1,lon1,lat2,lon2,r)
    real(8) :: spherical_distance
    real(8),intent(in) :: lat1,lon1,lat2,lon2
    real(8),intent(in) :: r
    
    !spherical_distance = r * acos( sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon1-lon2) )
    spherical_distance = r * acos(min(1.0d0, max(-1.0d0, sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2))))
  end function spherical_distance
  
end module nh_case6_mod
